Time structure of muonic showers 
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An analytical description of the time structure of the pulses induced by muons in air showers at 
ground level is deduced assuming the production distance distribution for the muons can be obtained 
elsewhere. The results of this description are compared against those obtained from simulated 
showers using AIRES. Major contributions to muon time delays are identified and a relation between 
the time structure and the depth distribution is unveiled. 
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Q ■ One of the highest priorities in Astroparticle Physics has become the understanding of the origin, composition and 
(-7- ' energy spectrum of Ultra High Energy Cosmic Rays. Cosmic rays have been observed for a long time with energies 
known to exceed 10^*^ eV but the remarkably low flux at these energies makes their detection difficult. The air 
' showers developing in the atmosphere from the interaction of these particles can be sampled at ground level with 
, arrays of particle detectors. Alternatively the fluorescence light emitted by the nitrogen in the atmosphere can be 
collected by optical systems . Both types of experiments have confirmed the existence of cosmic rays of such high 
] energies although there are discrepancies between the cosmic ray fluxes inferred by the different techniques 0, . 
• The study of extensive air showers of the highest energies is difficult because the atmosphere is part of the detector 
] and the interpretation of measurements is very indirect. The air shower is a complex chain of interactions resulting in 
fS| ■ a very large number of particles from which only sparse data samples are taken. From this sample, the properties of 
T— I ' the primary particle must be deduced. We must rely on simplified models, or on complex numerical simulations, to 
infer these properties. These simulations have intrinsic and unavoidable uncertainties arising from interaction models 
that require untested extrapolations of low energy accelerator data. One of the complications arises because it is 
[ often very difficult to discern physical effects due to shower fiuctuations or to primary composition from effects due 
■ to uncertainties in the interaction models. It is hoped that through redundancy in the observations, both the high 
energy models and the primary particle properties can be constrained. 
Q . Currently, the most important step in this direction is the development of the Auger Observatory. The southern 
?H ■ part of the Observatory is now under construction in Argentina and will combine the fiuorescence technique with a 
c/3 : particle array of water Cerenkov tanks covering a surface area of 3000 km^ . The signal from the tanks is recorded 
with Flash Analog to Digital Converters and relative timing is synchronized using GPS technology to 10-20 ns. The 
^ time structure of the shower front is of crucial importance because the shower arrival direction in an array of particle 
■""j ■ detectors is determined from the relative arrival times of the signals. This is relatively easy since the shower front 
rS ] can be approximated by a plane but clearly there are curvature corrections to the front and the time structure of the 
^ . signal and its fiuctuations will have an impact on the precision of the reconstruction procedure. Understanding the 
■ " " ' time structure of the signal is most important for controlling the angular resolution. Last but not least, preliminary 
observations at the Auger Observatory are already revealing that the time structure of the tank signals carries much 
information 6] and it is hoped that this new information can be used to answer fundamental questions such as primary 
composition. 

As most high energy cosmic ray detectors are, or have been, arrays of particle detectors, the time structure of the 
shower front has been studied before [^.Isllgl lTollTlLIT^ . However most work is either experimental or phenomenological 
and there is not a great deal of knowledge about the time structure of signals in air showers from the theoretical side 
or of its relation to the shower parameters. In this work we study the time structure of the muon signal in air showers. 
Our work is based on ideas that have been introduced for the description of muon densities in inclined showers 
and is a natural extension of these ideas. Similar ideas have been developed in [Hill [H. 

Muons are naturally quite energetic when they reach the ground: otherwise they would have decayed in flight. As 
they have long radiation lengths they do not suffer many interactions in their travel path from the point where they 
are created to ground level. As a result if muons are present in the signal they deflne the earlier part of the shower 
front, i.e. its onset which is usually chosen to deflne the arrival time to be used in the reconstruction of the direction. 
Also muons have a much flatter lateral distribution than photons or electrons and they always dominate the signal 
at ground level for sufficiently large distances from the shower axis. Moreover, for inclined showers, they dominate 
the signal everywhere at ground level because the electromagnetic part of the shower, generated by the decay of 
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neutral pions, is absorbed high in the atmosphere. Therefore, the description of the time structure of muons in air 
showers is of important practical interest for the reconstruction of air showers from measurements made in extensive 
air shower arrays. Understanding the muon part of the shower is not the complete story because there are important 
contributions to the shower signal due to electrons and photons. We will leave the study of the time structure of 
electrons and photons for future work. In this article we neglect the effect of the earth's magnetic field, although it 
is known that magnetic effects are important for the description of muon densities at high zenith angles. We have 
checked that the time structure of muons is not affected by the magnetic field up to zenith angles of 80 degrees and 
more, depending on the strength of the magnetic field. Although our results are general, and we expect that they 
will be valid at all zenith angles, we will concentrate on the large zenith angles range, where the muon signal is more 
important. 

The article is organized as follows. In Section II we present the main assumption on which this work is based and 
we develop an analytical description of muon energy distributions in air showers as a function of transverse distance. 
These will be used in the next section when the time structure of the shower front is addressed. In Section III we 
describe the two main sources of time delays associated with muons and we present a model for their description, 
discussing how it compares with the results obtained directly using simulations and the relevance of the magnetic 
field effects. In Section IV we apply the method to give general results about time distributions in air showers and in 
Section V we summarize and conclude. 



II. TOWARDS AN ANALYTICAL DESCRIPTION OF MUONS IN AIR SHOWERS 



A. Factorization of Muon Distributions 



Our starting point will be an approximate model for the description of the production of muons in air showers. As 
the shower develops, muons are in general produced at a range of altitudes with energies spanning a wide range and 
carrying some transverse momentum. We will assume that the muon distribution factorizes completely into a product 
of three independent distributions, the energy fi{Ei), the transverse momentum /2(pt), and the distance to ground 
level h{z) 
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where Ei and pt are the muon energy and transverse momentum at the production point, and z is the distance 
from this point to ground level measured along the shower axis. The relation between z and production altitude 
is straightforward but requires taking into account the curvature of the Earth for large zenith angle showers. The 
transverse momentum pt is defined with respect to the shower axis. Here the three distributions are normalized to 
1 and Nq corresponds to the total number of muons produced in the shower. After production we assume that the 
muons only undergo continuous energy loss and decay. This is known to be a reasonable approximation when the 
magnetic field is ignored. Hard catastrophic interactions, such as bremsstrahlung or pair production, are unlikely 
because the muon radiation length in air greatly exceeds the whole atmospheric depth even for the most inclined 
showers. The angular deviations of the muons reaching ground level due to multiple elastic scattering are negligible 
because the average muon energy at ground level is about 1 GeV for vertical showers and grows rapidly as the zenith 
angle increases to over 200 GeV for horizontal showers [isj . 

In Fig. ^ we show the geometry of the problem indicating the notation used. For a given distance to the shower 
axis, r, the actual distance traveled by the muon, is simply / = i/r^ + (z — A)^, where A = rtan^cosC. 6 is the 
shower zenith angle and C, is the polar angle measured in a plane perpendicular to shower axis. We choose it to be zero 
for the case in which / is a minimum. This system is natural for the description of the asymmetries of the problem 
and turns out to be most convenient. 

The distribution of production distance for the muons, h{z), depends on the details of the hadronic model, on the 
primary energy and on the primary chemical composition. We can parameterize it from Monte Carlo simulations. 
In this article we will use the results obtained with simulations run with the AIRES code 16] for an observation 
altitude of 1400 m, that of the southern Auger observatory, and using the QGSJET model for hadronic interactions 
|17| | . For zenith angles above 60° a Gaussian or a log-Gaussian distribution usually reproduces the production profile 
sufficiently well, but the distributions can be also described numerically if necessary. In Fig[21we show the distribution 
of muon production distance as obtained from AIRES and our fit to a log-Gaussian distribution. In Table we show 
the averages (logj^o zq) and widths (ctz) of log-Gaussian production distance fits performed for a number of showers at 
different zenith angles. 

Equation n assumes that the transverse momentum distribution is independent of both the parent pion energy and 
its production altitude. Quite generally, this is equivalent to assuming that the transverse momentum distribution 
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FIG. 1: Schematic representation of the geometry and notation used for obtaining the arrival time distribution of muons in 
air showers. To describe the position at ground level of the muon signal we will use a cylindrical coordinate system (z, r, ^) 
with the polar axis along the shower direction. Coordinates in the transverse plane r and C, are defined in such a way that 
^ = corresponds to the earliest part of the shower front. The dotted line represents the distance from first interaction (star) 
to the muon production point (dot). The dashed line corresponds to the path traveled by the muon, I, from which the time 
interval can be deduced. Since muons are observed as they reach ground level, observation at different angles (" will correspond 
to different time delays because the path to ground is shortened by A (which is measured along shower axis). 
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TABLE I: Results of fits to log-Gaussian functions (in decimal base) to the production distance distributions obtained with Aires 
for proton showers and an observation altitude of 1400 m for different zenith angles. The last two entries are the logarithmic 
average logjQ(2o/km) =< logj^Q(2/km) > and the width, az. The second entry is zo the distance associated to the logarithmic 
average. 

is a universal function which depends only on the hadronic interactions. When magnetic effects, muon interactions, 
and multiple elastic scattering are neglected, the transverse momentum distribution of the muons at production is 
the combination of that of the parent pion and that introduced through pion decay which is below 30 MeV/c and can 
be then ignored when compared to the typical momenta of hadronic interactions ^ 200 MeV/c. Fig. O illustrates the 
transverse momentum distribution of all the muons that reach ground level. The most distinctive feature that can 
be observed is a relatively sharp cutoff for large pt. Particularly for low zenith angles this is close to an exponential 
cutoff. A good approximation for f2{pt) is given by 

Mp,) ^ = Bp\ exp(-p,/Q), (2) 

iVo apt 

where A ~ 1, Q ~ 170 MeV/c is the characteristic momentum associated with hadronic interactions, and S is a 
normalization constant. 

The distribution obtained is actually similar to the transverse momentum distribution in hadronic interactions as 
could be expected if this was the only source of transverse momentum. The main corrections to this simplification 
arise through the cumulative effect of transverse momenta acquired in multiple scattering of parent hadrons before the 
final hadron (usually a pion) decays into a muon. This cascading effect contributes to enhance the high-pt tail of the 
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FIG. 2: Production distance for proton showers of 10^^ eV of 60° and 80° zenith angles as obtained using AIRES (histograms). 
Also shown is the result of a fit to a log-Gaussian distribution (continuous lines). The observation altitude was set to 1400 m. 



momentum distribution, but fortunately it has little effect on the shape of the bulk of the distribution. Indeed a small 
dependence of the average transverse momentum with production altitude can be understood, most likely because of 
this cumulative effect. As the width of the pt distribution greatly exceeds the magnitude of the shift in transverse 
momentum due to its correlation with production distance, one can effectively ignore this dependence which plays a 
secondary role. 

Finally the energy distribution of the muons is known to be most affected by muon decay which, depending on 
the pat h length traveled by the muons from the production point to ground level, affects muons of different energies 
|l3lll8| . One can, as a first approximation, use a fixed energy distribution at production, fi{Ei), which then evolves 
because of continuous energy loss and decay in flight, to give the distribution at ground level, which is experimentally 
accessible with an air shower array. The most important effect is that there is a strong suppression of the distribution 
for muon ground energies below a low energy cutoff which is controlled by the travel distance for the muons. Clearly 
this cutoff is heavily dependent on zenith angle because the average distance to production increases rapidly from 
1 km to 200 km as the zenith angle changes from vertical to horizontal. 

Eq. n ignores the transverse position of the parent particles (mostly pions) that decay into the muons. This is 
equivalent to assuming that these particles travel along the shower axis. A simple argument shows that indeed pions 
are confined to a relatively narrow cylinder compared to the transverse distances of interest in high energy air shower 
experiments at ground level. A pion of energy E and transverse momentum pt, will travel at an angle from the shower 
axis given by sm9 ~ cpt/E. The pion will on average travel a distance I = jcttt = E/ {mTrC^)cTT^ before decaying, which 
corresponds to a transverse distance — I sin 9 = r^^pt Ittit^ . Given that the pt distribution is suppressed for large 
Pt^ we can see that for A = 1 over 59% of the pions will have a transverse distance less than — r^2(3/m^ ~ 22 m 
because of decay, which can be neglected compared to the transverse size of an air shower at ground level. We also 
see that is independent of the pion energy. A simple calculation shows that taking into account the transverse 
distance increase due to multiple cascading implies less than a 25% correction to r.^- 20]. 

Eq. ^is made simple by ignoring all correlations between production altitude, energy, and transverse momentum, 
but it nevertheless has considerable predictive power. The fact that these correlations play a subdominant role in 
what it is to be described can be understood in terms of the relative importance of distribution widths with respect 
to correlation effects. Although some correlations are indeed present, we will show that the assumptions made, which 
ignore them, predict correctly all the major characteristics of the time structure of the muon signal in air showers 
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FIG. 3: Muon transverse momentum distribution at production for a proton shower of energy 10^^ eV and 60° zenith angle 
as obtained using AIRES. Only muons that reach an observation altitude of 1400 m are considered. The distribution of all 
muons (top curve) is compared to that obtained in different radial bins, [100,300] m, [30,100] m, [10,30] m (from second top to 
bottom). 



that we are going to describe. This can be considered sufBcient justification to ignore them. 



B. Energy distributions at different r 

For the subsequent discussion of arrival time distributions, we will need details about the energy distributions at 
fixed transverse position, r, corresponding to typical experimental situations. These distributions can be reproduced 
approximately by considering a simple energy spectrum at production and evolving the distributions accounting for 
energy loss and decay. We will assume fi{Ei) to be given by a power law 

ME,) ^ = AEpeiE, - mc^), (3) 

where mc^ is the muon rest energy and Q(x) is the Heaviside's step function. The parameter 7 can be fitted to Monte 
Carlo simulations. Results obtained in such way indicate that 7 has a very weak dependence on zenith angle. We 
have found 7 = 2.6±0.1 to be a reasonable value for all 9. 

We assume that the muon energy loss is constant dE/dx — —k where k — 2 MeV / (g cm^^) and take for simplicity 
a constant density atmosphere. In this case the muon energy can be simply expressed in terms of the path traveled, 
I 

EH) = E, - pkl, (4) 

where Ei is the muon energy at the production point. In addition the number of muons decreases with I due to muon 
decay. For a muon with initial energy Ei, losing energy in an uniform density atmosphere, it can be shown that the 
survival probability is 
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where the spectral slope k = mc^ / (crpk) ~ 0.8 for p = 1 x 10~^ g cm"^ ^M- instead of a single muon, we consider 
a muon energy spectrum at production, Eq. gives the energy spectrum after traveling a distance I in terms of the 
production spectrum. It can be seen that there is an effective low energy cutoff for muons that are produced with 
less energy than the corresponding loss during travel time. The low energy behavior of the ground energy spectrum, 
here characterized through k, is dependent on the muon energy loss and the atmospheric density. 

It is now straightforward to include energy loss and decay to obtain an approximate energy spectrum for the muons. 
We take the spectrum at production as given in Eq. Inland multiply it by the decay probability after traversing a path 
length I. The energy spectrum of the muons becomes 



dN 



anoE; 



Ei — pkl 



Q{E, 



pkl). 



Equivalently, the same distribution in terms of the final muon energy Ef is given by 

dN . ,„ . .V E 



dE 



f 



ANo{Ef+pkiy 



f 



Ef + pkl 



e{Ef~mc^). 



(6) 



(7) 



The energy spectrum discussed so far does not address the correlations between transverse distance and energy. 
The factorization assumption (Eq.^ allows us to explore the r-dependence of the energy spectra in some detail. This 
involves the transverse momentum distribution which is peaked at ~ Q. 

According to the assumptions discussed in the introduction, a muon produced at a distance z from ground level, 
with energy Ei and transverse momentum pt will arrive at ground at a transverse distance r given by 



Ptc 



^Ef - [mc^f 



cpt 
E,/ 



(8) 



where I = \/ [z — A)^ + is the total distance traveled by the muon and we assume yjE'f — {mc^Y — j given that 
Ei > nic^ + pkl. Eq. |S1 relates the transverse position of the muon to the three independent variables, Ei, pt and z 
and as a result all that remains is to perform the appropriate change of variables. 

In general we can convert from the pt distributions to r distributions simply by considering the Jacobian of the 
transformation 
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The last term in the above expression can be neglected provided that z 3> r, which is often applicable. 

If we approximate the pt distribution by the parameterization given in Eq. [5]thc spectrum at a given r is given by 
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In Fig. ^ we show the energy spectrum at constant transverse distance, r, given by the above equation compared 
to Monte Carlo results obtained with AIRES for a 80° shower. We have used the values of k and 7 obtained from 
theoretical considerations and we fit the value of pkl to the overall energy distribution at ground. The best fit value 
corresponds to an energy loss of pkl ~ 21.1 GeV. An alternatively approach could be to fit the parameters k and 7 in 
Eq. |3to the overall energy distribution from the Monte Carlo data and use them, as effective parameters to predict 
the energy distributions at fixed values of r. The results obtained are not critically dependent on the specific values of 
the parameters used nor on the form of the distributions used as inputs for the calculations. In fact it is well known 
that the energy loss parameter k is not constant but increases logarithmically with the muon energy. By taking it as 
a fitted parameter we include this energy dependence in an effective way. An analogous procedure can be applied to 
get the correlations between other distributions, for instance the pt distributions for different r which are shown in 
Fig. El 

For high zenith angles the production distance distribution is sharply peaked and fixing z to zq is a good approxi- 
mation. In Fig.^the results have been obtained in this way. For low zenith angles this approximation is not justified. 
We can instead consider the above equations to apply for a given z and integrate them over the z distribution. For 
example the energy distribution at fixed r becomes then 
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FIG. 4: Final muon energy spectrum at constant transverse distance, r. Histograms show Monte Carlo results obtained with 
AIRES for 10^® eV proton showers at 80° zenith angle compared to our prediction for the spectrum indicated by the continuous 
lines. The parameters used for these results were k = 0.8, 7 = 2.6, pkl = 21.1 GeV. 



For close to vertical showers this method of performing the z integral reproduces the results obtained with the 
simulations with much better accuracy. 

A close look at the distributions shown in Fig. 0] reveals that the energy distributions become broader as the 
transverse distance r decreases. This can be understood in terms of Eq. 1111 which predicts in general three distinct 
regions for the energy spectrum. 

1. For Ef < pkl the behavior is dominated by the term corresponding to the muon decay probability, i.e. by decay 
and energy loss. The spectrum is close to a simple power law dN/dEf ~ Ej. 

2. For Ef > pkl the spectrum reflects both the energy and transverse momentum distributions respectively through 
the indices 7 and A. The behavior of the energy spectrum is again close to a simple power law with a modified 
index which combines the two parameters, dN/dEf ^ EJ''^^^^ . 

3. For pt > Q the exponential suppression assumed in the transverse momentum distribution appears as a cutoff 
in the energy spectrum. It is simple to show that the suppression takes place for energies Ef > The cutoff 
energy grows as r decreases and this makes the distributions broader. This suppression affects the average 
muon energy at a given r and is responsible for much of the well-known correlation between average energy and 
transverse position for the muons. 

The three behaviors can be appreciated in the energy distributions for small r. For large r however the cutoff in 
transverse momentum arises for a final muon energy Ef < pkl. These distributions display a direct transition from 
the Ej behavior to the suppression at the large energy limit. The second condition {Ef > pkl) takes place at energies 
which are suppressed and thus has little influence on the bulk of the spectrum. We can estimate the approximate 

Qc 

distance Tc above which this happens, by combining the first and third conditions to obtain r > Vc ^ — - ^ 1000 m. 

pk 

Only for r below this scale the spectrum will display two slopes discussed above and a cutoff at higher energies. It is 
interesting to notice that is a transverse length scale which is approximately independent of zenith angle. It can 
be seen that this scale is also responsible for changes in quantities related to the muon energy in showers, for example 
the correlation between average muon energy and transverse distance. 

III. TIME DELAY ASSOCIATED WITH MUON PROPAGATION 

The time distribution of the muon signal at ground level is important. To study this we will consider the delays 
associated with different mechanisms in muon propagation. We measure the time delays with respect to a plane front 
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moving in the shower axis direction and traveling at the speed of hght. According to the simphfications discussed in 
section II, the muon arrival time will be determined by a combination of two effects. First there is a simple geometric 
effect (geometrical delay), muons arriving at distance r from the shower axis will take longer, because they travel a 
larger path than those arriving at r = 0. In addition these muons, being of finite energy, are further delayed with 
respect to the reference front moving at the speed of light because they travel more slowly (kinematic delay) . We can 
thus express the total time delay as a sum of these two delays 



(13) 



Here tg is the geometrical delay and is the kinematic delay. 

The distribution of arrival times taking into account both the geometrical and the kinematic delays can be now 
obtained by taking a convolution of the distributions associated with the two effects 



1 d^N 
Nr dtdr 



{t)®e{t) = / g{t - t')e{t')dt' , 



(14) 



where g{t) and e{t) are the arrival time distributions due to the geometrical and kinematical delays respectively, both 
normalized to 1 and iV^ is a shorthand notation for dN/dr which is needed for this normalization. We now study the 
two delays separately. 

A. Geometrical delay 

For high energy muons we can neglect, in a first approximation, the time delay associated with the finite muon 
energy. In that case we would simply convert the path length difference to a corresponding time delay. For a vertical 
shower we obtain 
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+ — z 
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This relation implies that the delay of muons arriving at a given point on the ground, at a fixed transverse distance 
r, is only a function of the production altitude. The arrival time distribution of the muon signal can be obtained as 
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1 d^N 
Nr dtdr 



dz 1 d^N 
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where hjz) , d"^ N / dEidr[z] (given in Ea. llO|l and dz/dt are functions of z and therefore implicit functions of t through 
Eq.E|[2ll. It can be shown that the z dependence of the integral of Eq.fTUIover Ei is l^''^ [1 — r^/P]. The first factor 
corresponds to an attenuation of the number of muons through decay 18] and besides the normalization has little 
importance for very inclined showers. Moreover since the distance distributions obtained with AIRES only account 
for the muons that reach ground level, it can be argued that the fit takes care of it. As a result the functional form 
of g{t) is given in implicit form through 



g{t) (X -Hz)^ 



r 



(17) 



where for large zenith angles the last factor can be dropped. 

It is relatively easy to understand that for inclined showers the geometrical effect implies that there is an asymmetry 
in the arrival times. The different path lengths traveled by the muons for different polar angles in the transverse plane 
(C) (see Fig. ^ induce different delays. We can take this into account by changing the relation between z and r of 
Eq. to include the excess path traveled by muons at different angles 



tgiz - A) = i (V(^-A)2 + r2 - (z - A)) 



(18) 



where A = rcosCtan^, and tg is measured with respect to the shower front plane. Now we still use Eg. 1161 to get 
the time distribution but dz/dt must be obtained from Eg 1181 Alternatively we can make the shift z ^ z + A in the 
argument of h(z) 



1 d'^N , , ^.dz 
'^'^-'N-rm-r°'-'^'^''^dt 



r2 



(19) 
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and still use Eq. ^]for dz/dt. As in our convention z is measured along the shower axis from ground to production 
point, a given dz/dt from Eg. 1151 must be associated with a production point which is shifted by A to ensure that 
distance to production corresponding to a point (r, C) is precisely z. The above equation introduces, in a simple way, 
important asymmetries in time distribution for constant r but different C, angles. 

In Fig. |5l we show an example of the pulse shape obtained using the result of Eq. ^] compared to the results 
obtained using Monte Carlo simulation carried out with AIRES. The agreement is reasonably good although there 
are significative differences. The predicted signal has a faster onset, is somewhat broader and has a shorter tail which 
is suggestive of a delay underestimate. We will show below that all these effects can be mostly attributed to the 
sub-luminal muon velocity. 




FIG. 5: Arrival time distribution of muons at r =2000 m for a 10^^ eV proton shower and incident zenith angle of 80° obtained 
with AIRES, for two values of = 0, 180°. The analytical result (continuous lines) with which it is compared takes into account 
only the geometrical delay. The observation altitude is 1400 m. 



Once we have an expression for the time distribution we can obtain all relevant statistical quantities associated 
with it. For example the average arrival time at a fixed transverse distance is simply given by 



< t{r) >- 



dt t g{t) ^ / dz h{z + A) tg{z). 



(20) 



For example, for r <C z, as is the case in inclined showers, we can expand tg given by Eq. 1151 to get a simple 
expression relating the time delay and the transverse distance 



< t{r) >= 



2c 



dz 



h{z + A). 



(21) 



For h{z) we can use a log-Gaussian distribution centered at zq with a standard deviation az, which is a good approx- 
imation for inclined showers (see Fig. ^ 



h{z) = — exp 
z 



1 

'^(logio(^/^o))' 



(22) 



where C is a normalization constant. When we substitute the above distribution into the expression for the average 
geometric time delay we get the following analytical expression 



< t(r) > = 



1 e 



2c ZQ- A 



(23) 
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where a = In(lO) (7^. The above expression gives the time curvature as obtained using the average arrival times of 
the signals as a function of the mean production depth. This curvature is directly related to the naive geometrical 
radius of curvature given by 1/2 r^/zg- We notice that there is a correction factor (involving the width of the muon 
production function distribution) between the naive geometrical curvature and the time curvature. 

We can also obtain the RMS value of the arrival time distribution of the muons at ground level as a function of r, 
zq and u. The result is 

o,(r) = V<t^>-<t>^ = ^ -^Vl-e-^ ~ ^- (24) 

ZC Zo — IC Zq — £\ 

This illustrates that the width of the arrival time distribution is proportional to the mean value, and hence it also 
rises as for a given shower. If is small the proportionality factor becomes simply a as indicated in the last 
equality of the above equation. 

In the same way we can obtained an estimation for the asymmetries in the geometric delay. We will define the 
asymmetry as 

A =< t{r, C = 180°) > - < t{r, C = 0°) > (25) 
Then a simple calculation using the distributions introduced before shows that 

r^tan^ ,„2 

A = 5— e^'" . (26) 

The asymmetry grows rapidly with r, as r^. Also, notice that at fixed r the average time can be approximated by a 
dipole type formula, < t >= a + bcos(. 




FIG. 6: Kinematical delay distributions for r ~ 10 m, 1 km, and 10 km obtained analytically. The parameters used are k = 0.8, 
A = 1, Q = 0.17 GeV/c, and 7 = 2.6. 



The main features of the arrival time distributions can be described with this simple model. As obtained in the 
above expressions, both the time delay and the width of the distribution grow as the transverse distance increases. 
The asymmetries obtained for different angles ( are in agreement with the simulation results. It is interesting to 
notice that the arrival times for small polar angles {( ~ 0°), corresponding to the part of the shower front that 
touches ground earliest, are further delayed with respect to the late arriving part {( ~ 180°). We can expect different 
curvature corrections for different polar angles. The asymmetries are largest for intermediate zenith angles, according 
to this simple model. In the actual experimental situation, typically involving an array of particle detectors, the 
curvature of the front is sampled at different angles ( so that systematic corrections to the arrival times and curvature 
of the front can be expected from this effect. 
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The results in Eq. 1241 may be useful when considering the weight to be given to a given time in the reconstruction 
of the arrival direction. These results can be significantly improved by considering the delay time due to the muon 
velocity which we now describe. 

A feature can be deduced from Eg. 1 191 which, in principle, opens up a new way to extract relevant information from 
the time structure of the shower front. It implies that a measurement of the arrival time of muons at a given distance 
to the core can give information on the production profile function of the muons. Important shower properties such as 
the shower maximum, which is related to the maximum in the muon production function, can in principle be obtained 
from it, to the extent that this approximation is valid. 



B. Kinematic delay 



There is an additional delay to be considered because muons with finite energy travel with a smaller velocity than 
that of light. This introduces the additional time delay te- For instance a 1 GeV muon traveling 40 km in vacuum 
would be delayed 750 ns with respect to the light front. In a medium such as air, muons have continuous energy loses 
and these delays are more complex to estimate. For this reason we have modeled the energy distributions in Section 
II-B, taking into account energy losses and muon decay. 

Neglecting multiple scattering and magnetic effects, muons travel a distance I = ^/z^+r^ from production to 
ground level, where z' = z — A. With these simplifications the extra travel time taken by a muon of velocity /3c with 
respect to a photon is simply given by 



1 

te = - 
c 



dl' 



1 



1 



(27) 



We can substitute into the above equation the muon energy as a function of travel distance I and then integrate. This 
can be done approximately using Eq. 01 to obtain 



I 



(28) 



where Ef = Ei — pkl is the final muon energy when it reaches the ground. An interesting feature is revealed by this 
equation. It leads to the conclusion that there is a maximum time delay which is independent of I and is approximately 
given by mc^/cpk ~ 1.7 fis. It corresponds to a muon that loses all its energy after traveling the distance I. This 
upper bound to the kinematic delay increases slightly if we consider an exponentially decaying density atmosphere. 
If the final energy is much larger than the muon mass, mc^ <C Ej we can expand Eg. 1281 to obtain 



1 {mc^y 

2 cpk 



1 



1 



Ei — pkl Ei 



(29) 



When we take into account the muon energy spectrum, for instance in terms of production energy Ei, the arrival 
time distribution is then given by 



1 (fN 
Nr dtfdr 



1 d^N dE, 
Nr dEidr dt. 



(30) 



At this stage, we make some simplifying assumptions to obtain simple expressions and discuss qualitative effects. 
If we expand Eq. 1291 for muon energies that are large in relation to the energy loss in the atmosphere we obtain the 
following relation between time delay and muon energy 



1 {mc^fl 

2 cE} ■ 



(31) 



If we now assume a fixed pt we can relate Ei to r to get 

1 (mc^)^ 



CPt 



(32) 



where we have expanded I as a function of z' and r for r ^ z' . This is quite interesting because we obtain a kinematic 
delay which has a very similar expression to that obtained for the geometric time delay in Eq. 1231 Very roughly 
one could expect an extra correction of order m? /p^ due to the kinematic delay. Assuming that the pt distribution 
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is sharply peaked at about 170 MeV/c this would result in ^ 25% more curvature. When the geometric delay is 
compared to simulation results, a systematic effect of this nature and of similar magnitude is indeed observed as 
remarked before. However the approximations involved to obtain this result are rather poor. 

Fortunately we can obtain more reliable time distributions due to the kinematic delay through Eq. I3UI We first 
invert Eq. to get 



2 



tc 



Here tc is a characteristic time that involves energy loss and the atmospheric density 



2\ 2 -, 
mc \ 1 



, , — ~ 500 ns 
pk / cl 



1 km 



I 



(33) 



(34) 



We now have the ingredients to calculate the kinematic time delay distributions. We have only to substitute the time 
derivative and the muon energy and transverse distance distributions (for instance as obtained in Section II) into 
Eq. 1301 The final result is cumbersome and we do not write it down explicitly. 

As discussed in Section II-B the energy distributions have different behaviors depending on the transverse distance 
r in relation to the critical radius rc- We can consider the two cases separately to better understand the effect of the 
energy distribution on the time structure of the pulses. For r ^ Vc we can simplify the energy distributions in Eg. 1111 
to a combination of a power law and an exponential. In that case the time distribution is simply 



e{t) cx t 



' exp 



1 {mc ) r 1 

2 pklcQ ct 



(35) 



For r Tc the time structure is more complex because we obtain two different behaviors in the two regimes described. 
This can be seen through Eq. 1331 in which the energy is regarded as a function of the delay time. The two regimes 
correspond to the times well below and well above tc- These regimes also correspond to the limits Ef ^ pkl and 
Ef <^ pkl respectively. In these two limits the expressions for the energy in terms of t can be simplified to 



• for t < tc Ei 



for t > Ef 



1 {mc^ Y ^ 
2 



ct 

1 {mc^f 

2 pkct 



When these relations are used together with the energy distributions to obtain the time distribution through Eq. 1301 
the two limiting behaviors obtained are 



eit) (X 



7-A-4 

t 2 exp 



1 mc r 1 
^c cQ ^t 



if t > 



if t < tc 



(36) 



Although we have made some effort to give analytical parameterizations where possible, making several approxi- 
mations, we stress here that the complete procedure, without these approximations, can be implemented. Note that 
Tc is below the distance separation between detectors in the Auger Observatory. 



IV. RESULTS 



We can combine the geometric and kinematic delays through the convolution given in Eq. ^| In Fig. |7| we show 
the average pulse time delay for the geometric and the kinematic delay times. We see that at small distances to the 
core the delay is dominated by the kinematic delay. This may seem surprising as near the core the characteristic 
muon energies are larger and one would expect a lower kinematic delay. However, near the core the spread on energy 
is larger (see Fig. and the time delay is dominated by low energy muons. 

In Fig.|S|we show our analytical result for the arrival time distribution compared to the Monte Carlo. We can see 
that the agreement is very good both at low and large times. This is illustrated further in Fig. where we show 
the time structure of muons arriving at transverse distances r ~ 400 m, 1 km, and 2.5 km compared to the results 
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FIG. 7: Average time delay for muons injected at 10 km with a 7 = 2.6 spectrum. The lower curve shows the geometric delay. 
The dashed points show the kinematic delay. The upper curve is the sum of both. 



obtained through direct simulation. These results have been obtained using k = 0.8 as suggested by theory, and 
A = 1, 7 = 2.6, Q — 0.17 GeV/c and pk— 0.2 GeV km^^. Our results are encouraging and we can state that the main 
contributions to the time structure of muons in air showers have been identified and understood. 

As an application of practical importance for direction reconstruction in extensive air showers we have calculated 
the average and the width of the arrival delay time distributions for different zenith angles as a function of distance 
to the shower axis, r. For near- vertical showers the production distance distribution can not be approximated by a 
delta function, as discussed at the end of Section II. The final kinematical delay distribution has been obtained by 
integrating the time delay corresponding to a given production distance z over the z distribution 

e(i;r)= J h{z)e{t; z,r)dz. (37) 

In Figs. ^1 ^2 El El we show the average and RMS values of the time distribution for times obtained with AIRES 
for 10^^ eV proton showers and our analytical calculation for 60° and 80°. It can be seen that the agreement is good 
specially the larger the zenith angle. The magnitude of the asymmetry is well reproduced as can be seen from the 
figure. For the RMS there are some discrepancies both at large and small distances. We believe that this is mainly an 
effect of fluctuations due to shower development and also due to the thinning procedure adapted in AIRES. In Fi g.lTOI 
we also show the parameterization given in Ref. based on a phenomenological formula suggested by Linslev |19j|. 

We have also calculated the asymmetries in the arrival time distributions. These can be approximated by a dipole 
correction to the average over ^ in the transverse plane. In Figs. PPil 1151 we give the magnitude of the asymmetry as 
a function of r for different zenith. The comparison between the analytical and the simulation results are given in a 
separate graph for characteristic distances. 

We consider that an important application of this method will be to provide a description of the time distributions 
of muons in a given detector. Since the distributions used are continuous, at some point the number of muons that 
reach a given detector will have to be selected, and a sampling effect will arise because the flnite number of muons 
that pass through a given detector. This sampling effect must be considered using the muon number distributions at 
ground level which have been discussed in detail in ref. [T^. 
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FIG. 8: Example of arrival time distribution from Monte Carlo (histograms) for a 10^^ eV proton shower and 80° zenith angle 
compared to our analytical calculation. Muons arriving in the range 3.2 < logjg r < 3.4 (r ~ 2000 m) are selected for two ( 
angles in the transverse plane. 



Up to now the magnetic field has been ignored but it will affect the time structure of the signal in two ways. 
Firstly muon delays are modified because particles do not travel in straight lines. Secondly significant modifications 
of the time delays calculated so far are to be expected because the simple geometrical relation between r, Ej , z and 
Pt is known not to hold. The effect of the magnetic field on this relation is the basis of the model in ref. [13] used 
to describe de density profiles of muons in inclined showers. There it is shown that the relative importance of the 
muon magnetic deflections with respect to deviations due to their pt is given by the parameter a — zB/pt. The 
condition of a > 1 corresponds io 9 > 80° at the Auger site. Clearly when a <C 1 one can expect our description 
to be valid. We have also tested the validity of our description for angles corresponding to a > 1 at the Auger site. 
Contrary to what could be expected, magnetic effects do not appear to modify the shape of the time distributions 
in an important way for zenith angles below 84°. Therefore the results obtained in this article without the magnetic 
field are phenomenologically interesting also in the region a 1 and for the Auger site the range of application is 
zenith angles below ^ 84°. 

The model described here on its own has also ignored shower fiuctuations. Nevertheless fluctuations of the arrival 
times of muons in an air shower can be related through the work presented here to fluctuations in the distributions 
discussed here. It is reasonable to assume that time fluctuations will be mostly due to the fluctuations in the 
longitudinal development of the shower and in the energy distribution of muons arriving at a given ground area. 

V. SUMMARY AND CONCLUSIONS 

We have obtained a procedure that can be used to describe the time structure of the muon signal in extensive 
air showers. This procedure is based on two simple mechanisms, time delays due to the extra path length which 
are discussed as geometrical delays and those due to the sub-luminal velocity of the muons. The global delay can 
be expressed as a simple convolution of these two delays. The description combines the physical mechanisms with 
phenomenological descriptions of the energy, transverse momentum and position of the muons that are produced in 
an air shower with several approximations. The results obtained agree sufficiently well with those obtained with 
direct simulations for most distances of interest and for all zenith angles to be of a great importance for practical 
applications. This agreement allows us to state with some confidence that these two mechanisms dominate the time 
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FIG. 9: Example of arrival time distribution from Monte Carlo (histograms) for a 10^^ eV proton shower and 60° zenith angle 
and = 180° compared to our analytical calculation. Muons arriving in the range 2.6 < logj^g ^ < 2.8, 3.0 < logjQ r < 3.2, 
3.4 < logj^Q r < 3.6 as marked. 
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FIG. 10: Average time delay as a function of transverse distance for two different positions in the transverse plane corresponding 
to ^ = 0° and 180°, as marked, for lO'^^ eV proton showers at zenith angle 60°. The histograms show the simulation results 
using AIRES and the lines show our analytical calculation. The lower (dotted)line shows the result of the phenomenological 
parameterization given in )12^| . 



structure of the muons in air showers. 



Acknowledgments 



LC and AAW acknowledge the hospitality of the Center for Cosmological Physics, University of Chicago, where 
this work was initiated. We also thank Maximo Ave and James Cronin for many discussions. This work was partly 
supported by the Xunta de Galicia (PGIDIT02 PXIC 20611PN), by MCYT (FPA 2001-3837 and FPA 2002-01161). 
RAV is supported by the "Ramon y Cajal" program. We thank CESGA, "Centro de Supercomputacion de Galicia" 
for computer resources. We thank the "Fondo social europeo" for support. AAW acknowledges continuing support 



16 



(0 

V 10' 
10^ 



10 



1 



10-^ 

1 1.5 2 2.5 3 3.5 4 

log,„(r/m) 
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FIG. 12: RMS values for the arrival time distributions of muons as a function of transverse distance for 10^^ eV proton showers 
at zenith angle 60° and for two different angles i^. The histograms show the simulation results using AIRES and the lines show 
our analytical calculation. 
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FIG. 14: Asymmetry in the average time delay which is plotted as a function of the angle in the transverse plane for 10^ 
eV proton showers at 60° and fixed r as marked. The histograms show the simulation results using AIRES and the lines show 
our analytical calculation. 
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FIG. 15: Same as Fig.[T31for zenith angle of 80°. 
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